%
% Tip of the Iceberg code : Welfare and trade barriers
%

% Parameters
clear;

workdir = '/Dropbox/RnD_KK/tradecost_replication';

% path(workdir,path);
addpath(genpath(strcat(workdir,'/ReStat_code')))


params.K = 2;
params.gdp = ones(params.K,1)*260029;
% t = 0;         % Alternative baseline
t = .13;

 ttar = 0;       % Baseline
% ttar = .047;   % (A) Counterfactual

params.tc = fliplr(eye(params.K))*t; 
params.ttar = fliplr(eye(params.K))*ttar; 
params.t = params.tc + params.ttar;

% Multiplicative barriers matrix, cols=to, rows=from
tau = 0;        
% tau = .3326;  % Alternative baseline
params.Tauc = ones(params.K) + fliplr(eye(params.K))*tau;

tautar = 0;        % Baseline
% tautar = .054496; % (M) Counterfactual
%tautar = .080124; % (M) eq tau baseline
%tautar = .039679; % (M) eq gamma=6

params.Tautar = ones(params.K) + fliplr(eye(params.K))*tautar;
params.Tau = params.Tauc.*params.Tautar;

% Baseline
params.sigma = 4;                  % Mean Broda Weinstein (median is 2.2)
params.gamma = 4;                  % Simonovska Waugh (was 4)

params.phi = 0;
params.P = 1;                      % Number of products
params.draws = 1e4;                % No. of random draws         

% Country size
params.L = params.gdp;    
params.mu = 1;                    % Consumption share for monopolistic good

% Random productivity draws. Productivity x is Pareto, F(x) = 1-x^(-gamma)
% Bounded from [1,+inf>, costs are <0,1]

rand('seed',1);
r = rand(params.draws,1);
params.x = (1-r).^(-1/params.gamma);

% Fixed costs
cxforeign = .533;       % CALIBRATE to share of firms exporting 
cxhome = 0; 
params.Cx = ones(params.K)*cxforeign-eye(params.K)*cxforeign+eye(params.K)*cxhome;   % (KxK)

% a constant
s = params.sigma;
params.delta1 = (s/params.mu)^(1/(s-1))*(s/(s-1));
params.potential_entrants = params.L*.15;   % CALIBRATE so that home cutoff=1

% Check what is t relative to median price
markup = s/(s-1);
%fobprice_median = markup./median(params.x).*(1+median(params.x).*params.t./(s*params.Tau));
%fobprice_median = median(markup.*(1./params.x));
%jj = (params.t./params.Tau)./fobprice_median;
%disp('t relative to median price'); 
%disp(mean(jj(1,2)));


%pause;

% ---------------------------------
% Check parameters before we move on
% ---------------------------------

if (~all(all(params.Tau.^(params.sigma-1).*params.Cx >= repmat(diag(params.Cx),1,params.K))))
    error('Error: Too low export costs');
end;
disp('Parameters loaded');

% What is the price index if x = mean(x)?
% Use the P0 as starting value.
p = markup.*(params.Tau./mean(params.x)+params.t);

for n=1:params.K  
  P0(n,1) = sum(params.potential_entrants.*p(:,n).^(1-s)).^(1/(1-s));
end

if (any(isnan(imo_solve(P0,params)))) error('Cannot solve price index. Problem ill conditioned.'); end;

% -----------------------------
% Solve equilibrium
% -----------------------------

func = @(p) imo_solve(p,params);
options = optimset('Display','iter','MaxIter',1000,'MaxFunEvals',20000);
tic; [P1,fval,exitflag,output] = fsolve(func,P0,options); toc;
if exitflag<=0 error('fsolve does not converge'); end;

% Find x_x and pi simulaneously
pi = solve_pi_xx(P1,params);
if (pi==0) warning('Profits negative, set to zero'); end;

Y = params.L*(1+pi);                     % y is total income
x_x = cutoff_imo_x(P1,Y,params);

disp('Results');
disp('Exitflag'); disp(exitflag);
disp('Price index'); disp(P1);
disp('Export cutoffs'); disp(x_x);
disp('L and Y'); disp([params.L Y]);
    
    % Find number of sellers and sales (simulated)
    for i=1:params.K           % Loop over source
        for n=1:params.K        % Loop over destination
    
           exporters = params.x>x_x(i,n);
           Entry_old(i,n) = sum(exporters);     % The # of firms exporting from i to n    
           price = markup.*(params.Tau(i,n)./params.x(exporters)+params.t(i,n));
           fobprice_median(i,n) = median(markup./params.x(exporters).*(1+params.x(exporters).*params.t(i,n)./(s*params.Tau(i,n))));
           p25price(i,n) = prctile(price,25);
           p50price(i,n) = prctile(price,50);
           p75price(i,n) = prctile(price,75);
           tlowerthantau(i,n) = sum(price./(price-t)<1.3326)/sum(exporters);           
           meanttotau(i,n) = mean(price./(price-t));
           medianttotau(i,n) = median(price./(price-t));
           
           s_in = params.mu*Y(n).*price.^((s-1)*(params.phi-1)).*P1(n).^(s-1);             
           share = sum(exporters)/length(params.x);                     % Pct. of firms above cutoff
           integral = share*mean(s_in);
           Sales(i,n) = params.potential_entrants(i)*integral;

           % Total labor allocated to producing goods in i exported to n       
           l_in = (s_in./price)./params.x(exporters);                  % Labor demand per firm 
           L_in(i,n) = params.potential_entrants(i)*share*mean(l_in);      

           %
           % Total t and tau tariff revenue          
           %          
           integral = share*mean(s_in./price);
           Qty(i,n) = params.potential_entrants(i)*integral; % Quantity sold
           Tot_t(i,n) = Qty(i,n).*params.ttar(i,n);
           Tot_tau(i,n) = Sales(i,n).*(params.Tautar(i,n)-1)./params.Tautar(i,n);           
        end   
    end

    Entry = repmat(params.potential_entrants,1,params.K).*x_x.^(-params.gamma);      % (KxK)    
    Avgsales = Sales./Entry;    
    
    % Check that total purchases equals total income
    disp('Total purchases and total income');
    disp([nansum(Sales)' params.L*(1+pi)]);    

    disp('Total manufacturing labor demand and supply');
    disp([nansum(L_in)' params.L]);
    
    disp('The number of products available in each country');
    disp(sum(Entry));
    
    disp('Imported products relative to domestic');
    disp(sum(Entry-eye(params.K,params.K).*Entry)./sum(Entry));
        
    disp('Entry matrix'); disp(Entry);
    disp('Sales matrix'); disp(Sales);
    
    disp('Share of firms exporting');
    disp(Entry(1,:)./Entry(1,1));
    disp('Cutoff matrix'); disp(x_x);    
    disp('Global profits');
    disp(nansum(nansum(Sales./params.sigma - Entry.*params.Cx))/sum(params.L));
    disp('Profit share - should be the same'); disp(pi);

    % Additive trade costs relative to price
    %idx = [1:17 19:22];
    idx = 2;
    t_no = params.tc(1,idx);
    tau_no = params.Tauc(1,idx);
    price_no = fobprice_median(1,idx);
    TC = (t_no./tau_no)./price_no;
    TCtar = params.ttar(1,idx)./price_no;
    disp('Trade costs relative to median export fob price'); disp(TC);
    disp('t tariff relative to median export fob price'); disp(TCtar);    
      
disp('Import share');
disp(sum(Sales-diag(diag(Sales)))./sum(Sales));

disp('Import share (imports/absorption) is 0.50 in data');
disp('Real wages'); disp(1./P1);

disp('Real wages incl tariff income'); disp((params.L+sum(Tot_t+Tot_tau)')./P1);

disp('Total t tariffs'); disp(Tot_t);
disp('Total tau tariffs'); disp(Tot_tau);

disp('Share firms with tau(t_A)<tau_M');
disp(tlowerthantau);

disp('Mean tau(t_A)');
disp(meanttotau);
disp('Median tau(t_A)');
disp(medianttotau);    


%%% ......................%%%
%%% TABLE  6              %%%
%%% ......................%%%

disp('  ');
disp('Table 6.............');


disp('Import share');
disp(sum(Sales-diag(diag(Sales)))./sum(Sales));

disp('Share of firms exporting');
disp(Entry(1,2)./Entry(1,1));